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We consider systems of many spatially distributed phase oscillators that interact with their neigh- 
bors. Each oscillator is allowed to have a different natural frequency, as well as a different response 
time to the signals it receives from other oscillators in its neighborhood. Using the ansatz of Ott 
and Antonsen (Ref. [23 ]1 and adopting a strategy similar to that employed in the recent work of 
Laing (Ref. [20]), we reduce the microscopic dynamics of these systems to a macroscopic partial- 
differential-equation description. Using this macroscopic formulation, we numerically find that finite 
oscillator response time leads to interesting spatio-temporal dynamical behaviors including prop- 
agating fronts, spots, target patterns, chimerae, spiral waves, etc., and we study interactions and 
evolutionary behaviors of these spatio-temporal patterns. 
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Many physical systems can be thought of as 
consisting of a large number of oscillating units 
that are distributed in space and coupled to 
neighboring units that are within some limited 
distance. The individual coupled units of such 
systems, moreover, can have non-negligible re- 
sponse times, and it is well known that delays 
can give rise to a set of possible behaviors that is 
significantly richer than would be the case with- 
out delays. Our work addresses two issues: (1) 
derivation of a macroscopic description for such 
systems, and (2) the possible characteristic be- 
haviors that may be revealed through study of 
such macroscopic descriptions. 



I. INTRODUCTION 

Systems of large coupled oscillator networks appear in 
many physical and engineering systems ll,]-[3j. Examples 
include synchronous flashing of fireflies [j] , pedestrian in- 
duced oscillations of the Millennium Bridge [E], cardiac 
pace-maker cells [6], alpha rhythms in the brain [7], gly- 
colytic oscillations in yeast populations jsj , cellular clocks 
governing circadian rhythm in mammals [9*], oscillatory 
chemical reactions [lQ|-[i2|, etc. 

Many previous studies of oscillator networks developed 
in the setting of network couplings on graphs with dif- 
ferent topological characterizations, such as small-world, 
Erdos-Renyi, and scale- free (e.g., Refs. [13]- [16]). Here 
we consider applications in which the oscillators are dis- 
tributed spatially, for example, when there is a row of 
trees each occupied with a large number of fireflies. In- 
deed, in the past decade studies of spatially distributed 
coupled oscillators have aroused much interest. An exam- 
ple is the chimera states (e.g., Refs. p/7j-[20]), in which 
there is a stable coexistence of both coherent and inco- 
herent states distributed in space. 

Another important aspect of the dynamics of oscilla- 



tor networks is that physical oscillators may have signifi- 
cant delays in their response to signals and these signals 
may also take a significant time to propagate. Studies 
of time-delay effects in the context of all-to-all coupled 
networks with a homogeneous distribution of time delays 
([21], [22]) show that interesting features such as bistable 
behaviors and multiple coherent states are induced in the 
presence of time delays. Reference |23], building on the 
machinery developed in Refs. [24]- [26], extends this line 
of study to a heterogeneous nodal response time distribu- 
tion. In addition, Ref. [27^ studies the dynamics of a one 
dimensional ring of spatially distributed and non-locally 
coupled oscillator network when the time delays are due 
to signal propagation between interacting oscillators. 

The problem studied in this paper is that of uncover- 
ing the spatio-temporal dynamics of a system of coupled 
oscillators with heterogeneous oscillator response times. 
We first give a microscopic description of the individual 
oscillators and their couplings. We then spatially coarse- 
grain this description and use the methods developed in 
Refs. [20| and [23| to derive a set of partial differential 
equations giving a macroscopic description of the system 
dynamics. Using our derived macroscopic equations, we 
then numerically explore the spatio-temporal dynamics 
and resulting pattern formation in both one- and two- 
dimensions. We find that a rich variety of behaviors are 
induced by the presence of time delay in the oscillator 
response. These include hysteresis, propagating fronts, 
spots, target patterns, chimerae, spiral waves, etc. 



II. FORMULATION 



We consider a system of N spatially distributed in- 
teracting phase oscillators with time delays between the 
response of an oscillator and the signal it receives. The 
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evolution equation of oscillator m is 

d ^ 

-^^rn{t) = CJm + ^ Kmn {sin[6>n(t - T^n) - ^m(^)]} 



f2 i^_i{e-^[^'»W-^"(*— )1 -c.c.}, 

(1) 



where Kmn is the interaction strength between oscillators 
m and n, which is assumed to be spatial in character (i.e., 
Krnn bccouies Small or zero if the distance between oscil- 
lator m and oscillator n is large), Tmn is the interaction 
time delay in the effect of oscillator n on oscillator m, 
and c.c. denotes complex conjugate. 

Assuming a separation in the scales of the macroscopic 
and microscopic system dynamics, we follow a path simi- 
lar to that employed by kinetic theory to reduce the study 
of a gas of many interacting molecules to a fluid descrip- 
tion. We begin by partitioning the continuous space into 
discrete regions centered at the discrete set of spatial 
points X, such that the domain of interest is Uxlx, and 
Ix r\ Ix' =0 for X ^ x' . The diameter of each region is 
|/^| ~ and the volume of each region is where d 
denotes the dimension of space. 

These regions are assumed to be small enough that 
Kmn ^ i^mi if oscillators n and / are in the same region 
Ix' ^ yet large enough that many oscillators ^ 1) 

are contained within each Ix' . Thus we can meaningfully 
define 



p(x') 
r{x',t) 



1 
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respectively as the local density and the local order pa- 
rameter in Ix' . In addition, for all m G Ix and n G Ix', 
we approximate ^ Kxx' • The summation in ([1]) can 
thus first be approximated as 
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(3) 

In all of what follows, we consider only the simple illus- 
trative case that Tmn = i-e., we suppose that the 
delay in the effect of oscillator n upon oscillator m is 
independent of n. This would, e.g., apply if the signal 
propagation time from n to m was very fast, but each os- 
cillator had a finite reaction time. Together with Eq.(|2j), 
Eq. dU can then be written as 

Y,w'K^^'p{x^)lm{e-^'-^'\{x\t- r^)}. (4) 

Since we assume Nj- ^ 1 for all x, it is appropri- 
ate to introduce a distribution function F{6^uj^x,r^t) 



proportional to the fraction of oscillators in Ix with 
e [0,0 ^ dO], uj e [uj.uj ^ dcj] and r G [r,r + dr] at 
time t. We furthermore pass to the limit of continuous 
space by replacing the discrete variable x by a new vari- 
able X which we now regard as continuous. In terms of 
this distribution, we introduce the marginal distribution 



g{uj,T, x) 



/ F{0,uj,r,x,t)dO. 
Jo 



(5) 



Here, note that since a;,r and x for any oscillator are 
assumed to be constant in time, the 6>— integral of F is 
time independent even though F itself depends on time. 
With Eq. (j5j), the quantity r in Eq. (j2j) becomes 



r I-oo C T, X, t)e^^dOdudr 

I-oo lo'' ^(0, uj, r, X, t)dOdudT 
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F{0,uj,T,x,t)e'^dOdujdT 



r(x, t) 



1 



p{x) Jo 



OO pOO 



oo ^0 



(6) 



The overall system dynamics can be studied in terms of 
the evolution equation for F(6>, a;, r, x, t), 

BP B 

— + - {F{co + lm[^{x, t - r)e-'0]}) = 0, (7) 



where 



r]{x,t) = J p{x')K{x,x')r{x' ,t)dx' 



(8) 



is Eq.(j4]) in the continuum limit, and the integration in 
([8]) is over the d-dimensional spatial domain. Referring 
back to our previous analogy to kinetic theory of a gas, we 
think of Eqs.([7j) and ([8j) as a kinetic description roughly 
analogous to the Boltzmann equation. 

To proceed we wish to reduce our kinetic description 
(O and (|8|) to a PDE (partial differential equation) sys- 
tem analogous to the fluid equations of gas dynamics. 
We do this using the recent work of Ott and Antonsen 
(Refs. j^-[^). We expand F in a Fourier series of the 
form 



Y,U{oJ,T,x,t)e'^' +C.C. 



I, l_n—i 

(9) 

As discussed and justified in Refs. [M] and [25,], we seek 
a solution in the form 



/n(a;, r, x, t) = d(cj, x, t - . 



(10) 
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Equations (j6j) to (j8]) then yield 
d 



dt 



a{uj^ x^t — r) -\- iuja{uj^ x^t — r) 

+ ^ [v{x,t - T)d?{uj,x,t - t) -7?*(x,t-r)] 

(11) 



(12) 



r{x,t) 



p{x) 



g{uj, T ^x)a*{uj^ x^t — r)dujdr ^ 



(13) 



where the star * denotes complex conjugate, and is 
written inside Eq. (p!3|) to emphasize its role as a dummy 
integration variable as compared with r's in the other 
equations. 

In what follows, we study an illustrative case corre- 
sponding to 



g{uj,T,x) = g{uj)h{r)po, 
Kix^x') = kqix — x')^ 



(14) 
(15) 



where g{uj)duj = h{r)dr = 1. Equation (p!4|) 
implies that the oscillator frequencies, locations, and de- 
lay distributions are uncorrelated, and that the oscilla- 
tor density po is uniform. Equation (fT5]) states that the 
strength of the coupling between oscillators at two points 
depends uniformly on their spatial separation. Further, 
in ([15]) we take q{x) to be suitably normalized, so that 
the constant k may be regarded as an overall coupling 
strength. With these assumptions, together with the 
transformation t ^ t -\- r in Eqs. (pTj) and ([T2]), and 
rewriting r' as r in Eq. (p!3|) . we obtain 



d_ 
dt 



a{uj^ X, t) + iuja{uj^ x, t) 
k 

+ - [7?(x,t)a^(cj,x,t) - 7^*(x,t)] =0, 
r]{x^t) = / poq{x — x')r{x\t)dx\ 



(16) 



(17) 



r{x^t) = / / g{(jj)a*{(jj^x^t — T)duj 

J lJ —oo 



h{r)dr. (18) 



In order to reveal generic expected behavior, we now 
further specify particular convenient choices for the fre- 
quency distribution, g{uj)^ the response time distribution, 
/i(r), and the spatial interaction kernel, q{x). 
We assume a Lorentzian form for g{uj), 



A/tt 



J-( !- 

27rz yuj — ujQ — 



iA UJ — ujo -\- iA 



(19) 



Assuming a is analytic in cj, we close the cj— integration 
path in (fT8|) with a large semi-circle of radius ^ oo in 



the lower half complex cj— plane. Thus we obtain from 
the pole of g{uj) dit uj = ujq — iA [see Eq. (p!9|) ]. 



ix^t) = la^ix^t-r)hir)dr 



(20) 

0, 

where Of ( x, t) = a{uJo — iA,x,t), and we have assumed 
(Ref. [2^) that, as lm{uj) — oo, a{uj^x^t) is sufficiently 
well-behaved that the contribution from the integration 
over the large semicircle approaches zero as ^ oo. 
Setting UJ = uo — iA in Eq. (p!6|) we obtain the following 
equation for the time evolution of a(x,t), 

d k 

—a{x,t)^{A^iuJo)a{x,t)^- [r]{x,t)a'^ {x,t) - rj* {x,t)] = 

(21) 

Our assumed form for the response time distribution 
h{r) is given by [23|, 



n-f3r,T 



(22) 



where An and /3n are defined by h{r)dr = 1 and 
rh{r)dr = T. Noting the convolution form of Eq. 
([2Q|) . we can re-express ([20]) as 



T d 



n^ldt 



n+l 



+ 1 r(x,t) = a*(x,t). (23) 



For the interaction kernel, we choose q{x) to be the 
solution to the problem. 



i2 



5ix). 



(24) 



For example, for an unbounded domain with boundary 
conditions q{x) ^0 as \x\ ^ oo, we obtain 



q{x) 



^exp(-^) ford=l, 
^i^o(^) for. = 2, ^^^^ 

_^exp(-M^ ford = 3, 



where Ko{\x\/L) is a zero order Bessel function of imag- 
inary argument. Using Eq. (|24|) . Eq. ([T7j) can be rewrit- 
ten by acting on it with the operator (V^ ~ z^)^ giving 



\/^r]{x,t) - j^r]{x,t) 



1 

"Z2 



por{x,t) (26) 



Thus we obtain a closed system of three PDE's in the 
independent variables x and t given by Eq. (|2T]) for 
a{x,t), Eq. ([23]) for r(x,t), and Eq. 1^ for r]{x,t). In 
the rest of this paper we study solutions of these equa- 
tions in one- and two- dimensional domains of size D 
with periodic boundary conditions. The parameters of 
this system are 

A,uJo,k,L,T,D, po,n. 
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By suitable normalization we can remove three of these 
parameters. We will do this by redefining r] and k to ab- 
sorb po and by normalizing time to and distance to 
L. This can also be viewed as using our original param- 
eter set with the choices A = l,L = l,po = l- In either 
case, our normalized PDE description becomes 



d 



—a{x, t) + (1 + iuJo)a{x, t) 
ot 



T d 
n^ldt 



n+l 



(27) 
(28) 
(29) 



In addition, in what follows we will only consider n = 
corresponding to h{r) = exp(— r/T). Thus, our 
reduced parameter set is 



uJo,k,T,D. 



(30) 



Before turning to the study of Eqs. (|27j)-(|29]), we 
briefiy comment on the analogy of the derivation of our 
evolution equations (|27|) - (|29|) to the derivation of the 
equations of gas dynamics from Boltzmann's equation. 
Substituting (p!Q|) into ([9]) and summing the geometric 
series (|a| < 1 is assumed for convergence), we obtain 



F(6>,a;, x, r, t) 



27r 



1 



+ -2|a|cos((9- V^) J ' 
_ (31) 

where a = |a| exp(— z?/^). It is shown in Refs. [2^ and [25| 
that, under very general conditions, the solution to our 
Eq. ([7j) relaxes to this form. In gas dynamics, the solu- 
tion to Boltzmann's equation, via the Chapman-Enskog 
expansion (Ref. [28]), is assumed to approximately re- 
lax to a local Maxwellian distribution whose velocity- 
space width is controlled by the temperature, and whose 
velocity-space maximum is located at the fluid velocity. 
In analogy with this situation, Eq. (|3T]) is peaked in 9 
(analogous to velocity space) at the location 6 = i/j (anal- 
ogous to the fluid velocity), and the width of this peak 
is controlled by |a| (analogous to temperature) with F 
becoming a delta function in 6> as |d| 1 (analogous 
to temperature ^0). In contrast to the derivation of 
gas dynamics from the Boltzmann equation, our relax- 
ation to (|3T]) is due to the phase mixing of many oscil- 
lators with different natural frequencies, whereas relax- 
ation to a local Maxwellian in gas dynamics is due to 
chaos in the collisional dynamics of interacting particles. 
Another difference is Uiat (|3T]) is an exact, rigorous result 
(as shown in Refs. [2^ and 0), while relaxation to a 
local Maxwellian in the derivation of gas dynamics is an 
asymptotic result in the ratio of the mean free path (and 
mean free time) to the macroscopic length (and time) 
scale. 



III. NUMERICAL STUDIES AND 
DISCUSSIONS 

A. ID propagating fronts, "bridge" and "hole" 
patterns 

The simplest solutions of our system, Eqs. (|27|) - (|29|) . 
are the homogeneous incoherent state solution (r = 
everywhere) and the homogeneous coherent state solu- 
tion (r = roe*^* where ro and Q are real constants). As 
shown in Refs. ^^-[23] for the case of globally coupled 
oscillators [corresponding to ^ in Eq. ([29]) ]. a 
distribution of interaction time-delays induces bistability 
and hysteretic behaviors. Figure [1] shows an example of 
the hysteresis loop in the \r\ — k plane for spatially ho- 
mogeneous states with cjq = 5,T = 1, which is obtained 
by solving Eqs. ([27|) to ([29|) with = r for the coherent 
solution r = roe*^^ 

We flrst consider a one-dimensional version of our sys- 
tem, Eqs. (|27|) - ([29|) . for a k value within the bistable 
region, /c = 12, and examine the evolution resulting from 
several initialized conflgurations with different spatial re- 
gions in the homogeneous incoherent and coherent state 
solutions. Results are shown in Fig. [2j Note that the 
flnal state is either coherent or incoherent depending on 
how large the initial incoherent region is. Thus, there ap- 
pears to be a critical initial size of the incoherent region 
beyond which the incoherent region takes over. Further- 
more, from Fig. [21 we see that the evolutionary process 
leading to this flnal state is by propagation of fronts sep- 
arating coherent and incoherent regions, and that these 
fronts propagate at an approximately constant speed. 
In addition to this initial example, we flnd a variety of 
other one-dimensional spatio-temporal behaviors to be 
reported in the following. 
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FIG. 1: Hysteresis loop for cjq = 5,T = 1. The upper and 
lower branches correspond to stable coherent and incoherent 
states. 



Next, we consider the dynamics as a function of the 
coupling strength k. Recall from Fig. [T] that there is 
a hysteretic region of coexisting coherent and incoher- 
ent states for the region k'^ < k < kc where /c^ = 10 
and kc = 14.5. Figure [3] shows the time evolution of 
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FIG. 2: |r(x,t)| for, (a) an initial configuration with a small 
part of the one- dimensional spatial domain in the incoherent 
state (blue) and a large part in the coherent state (orange), (b) 
a larger part of the spatial domain is initially in the incoherent 
state than that in (a), and (c) a still larger initial incoherent 
region, (uo = 5,T = 1, D = 100, k 12). 



\r{x^t)\ as a function of k. When the state is initialized 
with half (25 < x < 75) the domain in the homogeneous 
incoherent state and the remaining half in the homoge- 
neous coherent state, it is seen that if k is sufficiently 
close to /c^, then the incoherent region engulfs the co- 
herent region, while if k is sufficiently greater than /cc, 
the homogeneous coherent solution takes over, and by 
comparing Figs. 3(a) to |3(e)[ we find that the propaga- 
tion velocity decreases as k is increased toward kc. As 
k increases past k ^ 12, the simple propagating front 
phenomenon seen in Figs. [2] and 3(a)||3(c) is replaced by 
more complex behavior. For example, in Fig. |3(d)| we 
observe the formation of a "bridge" at A: = 13 (< /cc), 
i.e., a narrow stable coherent region sandwiched between 



two broad incoherent regions. This solution is apparently 
a long-time stable state. It develops as the two propa- 
gating fronts collapsing the coherent regions slow to a 
halt as they approach each other. We note further that 
the bridge has an amplitude which is smaller than that 
of the stable homogeneous solution, and the oscillation 
frequency is different as well (graphs not shown). Fur- 
ther, this bridge type solution persists for k > kc^ and 
can give rise to further intriguing dynamics like multiple 
bridges, as shown in Figs. 3(g)[ and even more vigor- 
ous behaviors of merging and re-creation of plateaus of 
coherent regions and bridges, as seen in Fig. [3(h)] Com- 
paring Fig. |3(f)| to [3(h)[ it is notable that a wide variety 
of evolutionary behaviors occurs within a relatively small 
range in /c, including the formation of single and multiple 
bridges, as well as collapse and re-creation of plateaus. 
Figure [3] studies the glassy-like behavior related to that 



seen in Fig. |3(h)| at a slightly different set of system pa- 
rameters. The figure shows plateaus of coherent regions 
(orange triangles in Figs. 4(a) and |4(b)[ ) and bridge-like 
patterns (yellow stripes), connected through dynamical 
creation, merging and re-creation of such structures un- 
til the system eventuall y evo lves into the homogeneous 
coherent state. Figure 4(c) shows the phase evolution 
inside the plateau region (orange triangle) of Fig. |4(a)^ 
centered at t ^ 420, x ^ 50. Figure [4(d)] shows the phase 
evolution corresponding to the four-bridge-structure be- 
tween the top of Fig. 4(a) and the bottom of Fig. |4(b)| 
(700 < t < 1300). We note that within a plateau, the 
whole region oscillates roughly homogeneously (see the 
nearly parallel evolving fronts in Fig. |4(c) ), and each 
bridge pattern functions as a sink of incoming waves (see 
the zig-zag-like pattern in Fig. |4(d)[ ). Further important 
dynamical characteristics during thi s vigo rou s glas sy-like 
transition state are revealed in Figs. |4(e) and |4(f)[ which 
show |r| and (where r = |r|e*^) respectively at t = 148. 
We see that there a re m ultiple hole-like patterns (deep 
dips in |r| in Fi g. |4(e)[ ), at which the phase changes 
sharply, (see Fig. |4(f)[ and note that the changes in phase 
for the outer two holes appear to be virtually discontinu- 
ous, as discussed in more detail shortly). In co mpari son, 
for the multiple-bridge region at t = 1200, Figs. 4(g) and 
|4(h)| show that both |r| and change smoothly in space. 

Figure [5] shows the dynamical characteristics associ- 
ated with the hole-like patterns in another setting where 
these patterns dominate and are not interspersed with 
other spatial features (like bridges and plateaus). The 
figure corresponds to the same parameters as those in 
Fig. |3(h)[ but initialized with different incoherent and 
coherent regions. Compared with Fig. |3(h)[ there is a 
relatively short time for the system to stay in the plateau- 
like regions, and instead of settl ing in the homogeneous 
coherent state solution as in Fig. |3(h)[ four distinct hole- 
like patterns emerge (black lines starting at t ~ 130 in 
Fig. 5(a)[ ). As time evolves, the two inner holes move to- 
ward each other and annihilate, while the outer two con- 
tinue to evolve, apparently becoming stationary. Note 
also that for the two merging holes, they approach each 
other at a faster speed when they are closer to each other. 
Examination of the phase evolution of the system (Figs. 
I and |5(e) ) suggests the center of each hole act as a 



source of plane waves, in contrast with the bridge solu- 
tion which acts as a sink (see Fig. |4(d)[ ). For the inner 
two moving holes, while each is characterized by a dip in 
magnitude (see Fig. |5(c) at t = 192), the dips decrease 
in magnitude as the two holes approach each other, with 
the relative phase difference on the two sides of the hole 
center close to being continuous (see Fig ]5(d)[ ). However, 
if the holes are stationary, e.g., the outer two holes in 
Figs. |5(c) and |5(f)[ each dip in |r| is close to zero, with 
the relative phase difference on the two sides being an 
essentially discontinuous slip of ±37r. A further observa- 
tion in the case of two stationary holes is that there is a 
bump in |r| half-way between them corresponding to the 
location at which incoming waves emitted from the holes 
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converge (see x = 50 in Fig. |5(f)[ ). 

In fact, when k ^ kc, the hole-hke pattern is a feature 
that shows up readily when two plane waves with a rel- 
ative phase difference of ±7r (or odd-multiples of them) 
collide. An example is studied in Fig. [6] where two waves 
of relative phase difference tt collide giving rise to a hole 
pattern. This observation is consistent with the relative 
phase differences observed at the two outer holes stud- 
ied in Fig. |5(g)[ Furthermore, although the hole pattern 
seems to arise only under relatively specific conditions, 
it is found to be pretty stable with respect to changes 
in parameters or small perturbations once it is formed. 
Finally, as shown in Fig. [71 we note that the hole core 
occupies a finite width and so is not a point singularity 
when T 7^ 0. This will be shown to have a close corre- 
spondence with the spiral wave in our two-dimensional 
study (sec. HITC]) . 

It is further interesting to note some similarity between 
our observations in the region k ^ kc and the intermit- 
tency regime of the complex Ginzburg-Landau equation 
(CGLE) (See for example, section III of Ref. [29], and 
section 2.5 of Ref. [30]). There, the CGLE displays sim- 
ilar glassy-like transition patterns characterized by large 
plateaus of coherent regions with hole-like patterns being 
continuously created and destroyed. However, there are 
also differences between the two systems. For example, 
the CGLE does not seem to have a close counterpart to 
the bridge pattern observed in our system, while more 
intricate dynamics of hole creation and destruction lead- 
ing to zigzagging holes and defect chaos have not been 
observed in our study. 



B. 2D propagating fronts and "bridge" patterns 

Figures [8] and [9] show the d = 2 counterparts to the 
d = 1 propagating fronts and associated features. Similar 
to what was previously done for d = 1, half of the system 
is initialized in the homogeneous incoherent state and 
the remaining half in the homogeneous coherent state, 
and they are divided by a sinusoidally wiggling bound- 
ary (Figs. 8 (a) I and 9(a)). Analogous to the d = 1 case, 
for d = 2, the homogeneous incoherent state and homo- 
geneous coherent state take over when k is sufficiently 
small or large compared to kc respectively. The most in- 
teresting behaviors again take place when k ^ kc- With 
k = 14.4 < kc, Fig. [8] shows the development of a sta- 
ble bridge solution. In contrast, with k = 14.8 > /Cc, 
Fig. [9] shows a surprisingly rich spatio-temporal pattern 
evolution. As in Fig. [8l the originally coherent half ap- 
parently starts to shrink into a bridge (see Fig. |9(b)[ ); 
however, as time progresses further, we see that coher- 
ent regions arise out of the originally incoherent regions 
to form new features (see also Figs. 3(g)| and |3(h)| in 
the d = 1 case), and these new features interact in a 
nontrivial two dimensional manner. For example, when 
two neighboring coherent regions get close to each other, 
they can form bonds and merge into each other: see the 



connections formed between bridge-like structures from 
t = 83 to t = 98; also see the coherent spot formed at 
the upper left hand corner at t = 245 and see how it 
merges into the bridge on its right as time progresses to 
t = 400. We also observe that, during the process of 
merger, bridge-like structures may also temporarily sep- 
arate and then re-connect: see the connecting bridge at 
the bottom right hand corner from t = 138 to t = 170. 
A further notable feature is the coherent spot on the top 
left hand corner at t = 561 (a target pattern in the phase 
plot as shown in the next section), which survives from 
t = 561 to the end of the numerical run. In the above 
reported numerical experiments we observe that both in- 
coherent and coherent regions coexist for a long time. 
We do not know, however, whether a coherent or inco- 
herent state ultimately will take over the whole domain 
at longer time. 



C. 2D Spots, spiral waves and target patterns 

Figure [To] shows the time evolution of both |r(x, t)| 
and sin[6>(x, t)] [ where r(x,t) = |r(x, t)| exp[i^(x, t)], and 
X = (x, y) in 2D ] when the system is initialized with a 
small random initial condition at each grid point, and 
the coupling strength is k = 15.5 (> kc = 14.5). As 
expected from our previous studies, when k > kc, coher- 
ent regions (|r| ~ 1) emerge from the initial incoherent 
state. Further, the phase plots show some distinct target- 
like patterns of nested closed surfaces of constant phase 
(see t = 40 and t = 217). As time progresses, coherent 
regions (red in the |r| plots) become dominant and only 
small islands of incoherent regions remain (blue in the |r| 
plots). Similar to our previous observation of propagat- 
ing fronts when k > kc (compare Figs. |9(g)| and |9(h)[ ), 
coherent regions can form in an originally incoherent re- 
gion (|r| <C 1). For example, see the figures from t = 139 
to t = 161, and especially from t = 195 to t = 225, where 
we see coherent regions (red/yellow) appearing and grow- 
ing in the interior of incoherent (blue) blob, eventually 
destroying it. As can be inferred by comparing the \r\ and 
sm{6) plots, small blue, dot-like features in the |r| plots 
represent phase defects in the complex amplitude (i.e., 
counter clockwise encirclement of such a feature leads to 
a phase change of either +27r or — 27r), and these blue dot 
features are commonly seen as spiral wave type patterns 
in the phase plots. When, as in the previously noted 
plots from t = 195 to 225, coherent regions take over 
from an incoherent patch, we also note that a number of 
phase defects result (which must be formed in opposite- 
spiral-parity pairs due to the conservation of topological 
charge); see t = 250. The isolated phase defects sub- 
sequently wander about, and some of them are seen to 
annihilate with others of opposite parity (see the two de- 
fects closest to the bottom of the picture at t = 267 and 
their evolution up to t = 293), or sometimes they are 
absorbed into an incoherent region (e.g., compare the \r\ 
plots at t = 195 and t = 217). Lastly, regarding the 
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speed of motion of spiral patte rns, we note that similar 
to the observation in Fig. |5(a)[ when oppositely charged 
spirals get close enough to each other, their speed of ap- 
proach becomes distinctively faster till they annihilate 
each other. 

In studies of the CGLE, the hole pattern and spi- 
ral wave pattern are analogous phenomena occurring in 
d = 1 and (i = 2, respectively. Indeed, the hole pat- 
tern and spiral wave pattern exhibit similar characteris- 
tics in our study. Both features are stable with respect 
to small changes in parameters, and exhibit similar dy- 
namical characteristics of approach and annihilation as 
described above. In addition. Fig. [TT] shows, in parallel 
with Fig. [71 that the central core of the spiral wave pat- 
tern occupies a finite area when T ^ 0. This is similar 
to the chimera- centered spirals noted in Refs. .32j . 



D. 2D pulsating pattern 

Another class of local coherent structures supported in 
the d = 2 case is shown in Figs. [T2land [T3l which shows 
a localized pulsating spot in an incoherent background. 
It is interesting to notice that oscillations of the mag- 
nitude and phase (which show up in the form of target 
patterns) of r(x, t) are not the same, with that of the 
phase oscillation being more irregular and more than an 
order of mag nitu de fast er than the amplitude oscillation 
(Figs. |12(d)| and |13(g) ). It is interesting to note that for 
the CGLE, stable pulsating patterns come only with the 
addition of a quintic term (see Ref. [33], and the later 
work Ref. [s^ and references therein). 



IV. SUMMARY AND CONCLUDING 
REMARKS 



rived a macroscopic PDF description for this situation 
[Eqs. (|27|) - (|29|) ]. The resulting macroscopic dynamics 
are found to exhibit a wide variety of pattern forma- 
tion behaviors. We characterized the possible behaviors 
roughly according to the hysteresis loop corresponding 
to bistable homogeneous incoherent and homogeneous 
coherent state solutions. Numerical studies show that 
the system behaviors for k sufficiently below/above the 
bistable /c-range are simple in that the homogeneous in- 
coherent/coherent state eventually takes over the entire 
domain. In contrast, for k in or near the bistable range 
the system can exhibit a variety of interesting spatio- 
temporal phenomena. These include propagating fronts, 
bridge patterns, hole patterns {d = 1), spiral waves 
{d = 2), spots, target patterns, pulsating patterns, etc. 

Finally, it is interesting to consider the role of time 
delay in contributing to the features that we observed. 
If there is no time delay (i.e., T = 0), there is no ho- 
mogeneous bistable behavior as observed in Fig. [U and 
the transition from the homogeneous incoherent state to 
the homogeneous coherent state is supercritical and takes 
place at kc = 2A. In this case, many of the interest- 
ing spatio-temporal phenomena that we have found for 
T > are absent. For example, when T = 0, the intri- 
cate ID glassy state transitions were not observed, and 
the system typically evolves relatively rapidly into either 
homogeneous incoherent or homogeneous coherent state 
solutions. The 2D waves arisen from topological defects 
are still present; however, for T = the system will be 
similar to the case of zero nonlinear dispersion in Ref. 
[l^, where the incoherent core remains a point defect 
but not a finite area as observed when T ^ 0. Thus finite 
response time introduces additional dynamics, leading to 
the large variety of behaviors observed. 



In this paper, we have studied the spatio-temporal dy- 
namics of spatially coupled oscillator systems where the 
oscillators have a heterogeneous distribution of response 
times. Using the results of Refs. [2^]-[26j, we have de- 
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FIG. 3: A comparison of the time evolutions of |r(x,t)| for different values of k where r is initialized with half of the interval 
at th e coherent state (25 < x < 75) and half at the incoherent state. Notice the difference in time scales of Fig. 3(g) and Fig. 
3(h) from other figures (cjq = 5,T = 1, D = 100; periodic boundary conditions are imposed). 




FIG. 4: (a,b) Glassy state of transition, formation of plateaus of coherent regions and hole patterns, and final evolution into the 
homogeneous coherent state. (c,d) Phase evolution in the plateau and multiple-bridge regions. (e,f) |r(x, 148) | and 0{x, 148). 
(g,h) |r(x, 1200)1 and 0{x, 1200) (cjo = 4, T = 1, D = 100, k = 10.3 (> kc = 10); initial condition: r is given by the homogeneous 
coherent state solutions for 25 < x < 75, and r = otherwise; periodic boundary conditions are imposed). 
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{e)e{x, t) (f)r(x, 1200) {g)0{x, 1200) 



FIG. 5: Formation and dynamical evolution of hole patterns, (a) |r(x, t)|. (b-d) Close-up views of four hole patterns with two 
inner traveling holes, (e-g) Close-up views of two stationary hole patterns, (cjq = 5, T = 1, D = 100, k = 14.8; initial condition: 
r is given by the homogeneous coherent state solutions for < x < 41 and 59 < x < 100, and r = otherwise; periodic 
boundary conditions are imposed). 




(a)\r(x,t)\ {h)0(x,t) (c)|r(a;, 200)| (d)6'(x, 200) 



FIG. 6: An example of the hole solution by collision of two plane wave solutions. The two waves meet at x = 50 with a tt 
phase difference, {uq = 5,T = l,i^ = 100, /c = 14.8 and periodic boundary conditions). The initial condition corresponds 
to a discontinuous r given by a right traveling plane wave solution with m = 3 (where the wave number is 2rmT / D) for 
< X < 50 and a left traveling plane wave solution with m = 4 for 50 < x < 100. Correspondingly, we observe from (d) that 
[6'(0,200) - 6'(100,200)] = 2n. 




X X 



(a)|r(x)| (h)e(x) 



FIG. 7: Finite width of (one-dimensional) hole core (cjo = 5,T = 1, D = 33.3, k — 14.8). 
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FIG. 8: Time evolution of |r(x,t)| of a d = 2 configuration initialized with half of the region at the incoherent state and half 
at the coherent state divided by a wiggled boundary with k = 14.4 (< kc = 14.5) (cjq = 5,T = 1,D = 100; periodic boundary 
conditions are imposed). 
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FIG. 9: Time evolution of |r(x,t)| from an initial configuration (a) with half of the region at the incoherent state and half at 
the coherent state divided by a wiggled boundary with k = 14.8 (> kc = 14.5). A comparison with Fig. [8]shows a much richer 
spatio-temporal dynamical pattern (cjo = 5, T = 1, D = 100; periodic boundary conditions are imposed). 
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FIG. 10: Time evolution of |r(x,t)| and sin[^(x,t)] (where r(x,t) = |r(x, t)| exp['i^(x, t)]) from random initial condition, (a) 
and (b) (cjq = 5,T = 1, = 100, k = 15.5 (> kc — 14.5); periodic boundary conditions are imposed). 
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FIG. 11: Finite area of (two-dimensional) spiral cores (cjq = 5,T = 1,D = 20, /c = 15). 




FIG. 12: Pulsating pattern: amplitude variation. Figures (a) to (c) show approximately one "period" of oscillation in |r|. 
Figure (d) shows the time variation of |r| at the center of the pulse; compare with Fig. [13] for oscillations in phase (cjq = 5, T = 
1, D = 100, k = 14.52; periodic boundary conditions are imposed). 
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FIG. 13: Pulsating pattern: phase variation. Figures (a) to (f) show the rapid time variations of the phase. Figure (g) shows 
the time variation of sin{0) at the center of the pulse (Parameters are as indicated in Fig. I12p . 



